Dynamical heterogeneity in soft particle suspensions under shear 
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We present experimental measurements of dynamical heterogeneities in a dense system of microgel 
spheres, sheared at different rates and at different packing fractions in a microfluidic channel, and 
visualized with high speed digital video microscopy. A four-point dynamic susceptibility is deduced 
from video correlations, and is found to exhibit a peak that grows in height and shifts to longer 
times as the jamming transition is approached from two different directions. In particular, the time 
for particle-size root-mean square relative displacements is found to scale as r* ^ (7A(/)^)~^ where 7 
is the strain rate and A(/) = — (j)c\ is the distance from the random close packing volume fraction. 
The typical number of particles in a dynamical heterogeneity is deduced from the susceptibility peak 
height and found to scale as n* (7A(/)^)~°'^. Exponent uncertainties are less than ten percent. 
We emphasize that the same power-law behavior is found at packing fractions above and below 
(pc Thus, our results considerably extend a previous observation of n* ^ 7"°-^ for granular heap 
flow at fixed packing below (j)c. Furthermore, the implied result n* ^ (t*)°"^ compares well with 
expectation from mode-coupling theory and with prior observations for driven granular systems. 

PACS numbers: 64.70.pv, 83.80.Kn, 05.20.Jj 



I. INTRODUCTION 

Disordered materials of all kinds are considered to be 
"jammed" if the relaxation time grows longer than the 
observation window, so that that the constituent par- 
ticles appear locked into a fixed configuration of nearest 
neighbors [TH3]. For example, supercooled liquids can be- 
come jammed by lowering the temperature; hard sphere 
colloidal particles can become jammed by increasing the 
density; macroscopic glass beads can become jammed by 
lowering a driving force below some threshold. No mat- 
ter what the material or set of control parameters, as 
jamming is approached it has long been assumed that 
the growing relaxation time is accompanied by increas- 
ing co-operativity in particle motion [4 . The closer to 
jamming, the larger the number of neighbors that must 
cooperate in order to rearrange and the less frequently 
this happens. 

It is now widely accepted that the rearrangement dy- 
namics are not continuous near jamming, but rather are 
spatially and temporally heterogeneous [5 . Intermittent 
string-like swirls of rearranging particles come and go in 
a background of less mobile particles. The four-point dy- 
namical susceptibility X4.{^) is a powerful tool for charac- 
terizing such dynamical heterogeneities [6l [7] . This func- 
tion exhibits a peak at a characteristic relaxation time, 
r*, and the peak height x% can be related by a counting 
argument to the number n* of particles in the fast-moving 
rearranging regions [8]. One of the central questions to- 
day, then, is the quantitative relationship between the 
respective growth of r* and of xX on approach to jam- 
ming. Expectations for various models are reviewed in 
Ref. [7]. For example a logarithmic connection is ex- 
pected for "collectively-rearranging region" scenarios. A 
power-law connection xl {r*)^ is predicted by mode- 
coupling theory, where A is the reciprocal of the mode- 



coupling exponent, 7; Ref. [3 particularly notes the val- 
ues A = 0.37 [9 and A = 0.40 [10]. A power-law connec- 
tion with A = 1 is expected for freely-diffusing defects. 
And more recently a value A = 1/2 was reported for a 
kinetically constrained model jamming model pT] , 

For colloidal hard spheres this issue was recently ex- 
plored in Ref. p^, which improves upon pioneering ob- 
servations [T3HI6] by covering an unprecedented density 
range near jamming such that the structural relaxation 
time increased by seven orders of magnitude. The data 
show that r* grows faster than a power law, and n* grows 
slower than a power law, in l/((/)c — (/>) as ^ approaches 
(j)c from below. The critical packing fraction (pc is close 
to, but possibly distinct from, random close packing. Ir- 
respective of the value, the conclusion is that n* grows 
logarithmically with r*. For other colloidal systems it is 
not yet known whether this relationship depends on the 
nature of the particle interactions, or whether it changes 
when the control parameter is temperature or driving 
rather than just density. In this paper, we report on 
dynamical heterogeneities for dense suspensions of soft 
Hertzian colloidal particles. In particular we measure 
the time and size scales, and determine how they grow as 
jamming is approached both by bringing the density to- 
ward (j)ci from either side, and also by lowering the strain 
rate. In addition to establishing the dependence of n* 
and r* on these two control parameters, we also show 
that the size and time scales are related by a power law. 
This result contrasts with Ref. [12 , but compares well 
with observations for macroscopic hard spherical grains, 
where the control parameters are density and fluidizing 
air speed p/7j or strain rate and depth into a flowing 
heap [H]. Here, as in Ref. ^ and Ref. [^, the dy- 
namic range in relaxation time is more than seven orders 
of magnitude. 




FIG. 1: (Color online) Properties of NIPA microgels as a func- 
tion of packing fraction, as controlled by temperature: (left) 
the elastic modulus measured by centrifugal compression, and 
(right) the particle diameter measured by dynamic light scat- 
tering [19]. 



II. EXPERIMENTAL DETAILS 

The system we study is a dense aqueous suspension of 
thermoresponsive N-isopropylacrylamide (NIPA) micro- 
gel beads [2Ql |2T] , synthesized with the Yodh group at 
Penn [22H26] . Dynamical heterogeneities for unsheared 
suspensions of such particles have been reported previ- 
ously in Refs. [24l [27l |28] , both below random close pack- 
ing as well as above - where aging effects are important. 
Here two different size particles are used, primarily about 
1 //m but also about 0.6 jam in diameter; for the former, 
the number density is 0.455//im^ and the viscosity of the 
suspending water is r]o = 0.01 g/(cm-s). Fig. [l] shows 
the diameter and Young elastic modulus E for the larger 
particles, obtained previously from dynamic light scatter- 
ing and centrifugal compression [19 respectively. Note 
that decreasing the temperature causes the particles to 
swell with water and become softer. The applied pres- 
sure needed to squeeze water from the gel is very large 
compared to the elastic modulus [19]; therefore, the par- 
ticles deform without deswelling and can be compressed 
to a known volume fraction above random close pack- 
ing 0c = 0.635, simply by lowering temperature. 

Previously we studied the shear rheology of these sus- 
pensions by a custom microfluidic technique, in which 
the velocity profile is measured at the mid-height of a 
tall channel for various packing fractions and for various 
pressure-controlled flows [29]. The channel is 25 /im wide, 
100 /im tall, and L = 2 cm long, fabricated of PDMS by 
soft lithography and bonded to a glass microscope slide. 
The suspension is forced through the channel using pres- 
surized air and inlet/outlet tubing of sufficient diameter 
that the imposed pressure drop AP occurs only along 
the length L of the channel and can be related to the 
local shear stress in the suspension as (j{y) = /S.Py/L. 
The local shear strain rate is found by numerical dif- 
ferentiation of the velocity profile, 7(7/) = dvx{y) /dy. 
For this, we collect video data with a Phantom CMOS 
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FIG. 2: (Color online) (a) An image taken from video data 
of particles, (b) The velocity profile across the channel, ex- 
tracted using PIV, and (c) the strain rate as a function of 
position. 



camera (1-10,000 fps) connected to a Zeiss Axiovert 200 
microscope with 100 x objective. An objective-cooling 
collar (Bioptechs) and cooling plate above the sample 
are controlled to about 0.1 C in order to vary the parti- 
cle volume fraction. An example video frame in Fig. [2^ 
displays bead-scale intensity variations, so that Particle 
Image Velocimetry (PIV) may be implemented with cus- 
tom Lab VIEW code. In practice, we break the image 
into strips 7 pixels wide and determine the velocity and 
strain rate in these regions. Example velocity and strain 
rate profile data are shown in Figs. [2]3-c. In this figure 
the packing fraction is = 0.643, which is slightly above 
random close packing (j)c = 0.635; therefore, the flow is 
somewhat plug-like and exhibits wall slip. We found that 
the resulting stress vs strain rate shear rheology could be 
collapsed onto two branches, by Olsson-Teitel [30] scaling 
with powers of the distance = \(j) — (j)c\ to jamming; 
the resulting exponents can be understood in terms of 
particle interactions [31] . 

The full data set from Ref. [29] consists of video plus 
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FIG. 3: (Color online) (a) Overlap order parameters and (b) 
corresponding dynamic susceptibilities, plotted vs delay time 
for several strain rates as labelled. The packing fraction is 
= 0.643, as in Fig. [2] where strain rates are seen to be lower 
towards the center of the channel. In (a) the light dashed 
curves represent Qt{r) for a selection of different start times 
t, and the heavy solid curves represent the average Q(t) = 
{Qt{r)) over all t. 
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FIG. 4: (Color online) Averages of the overlap order param- 
eter at time delay r*, plotted versus volume fraction, where 
each data point represents a different strain rate: (a) Qi and 
(c) Qo are the averages for the slow and fast regions, respec- 
tively, whereas (b) Q* = Q{'^*) is the average for the whole 
sample. None of these quantities is found to depend on vol- 
ume fraction or strain rate; their averages are indicated by 
the dashed horizontal lines with accompanying values. 



stress and strain rate profiles for packing fractions varied 
discretely between 0.5 and 0.7, and for strain rates varied 
continuously from 10~^ to 100 s~^. Since the systems are 
always under shear, steady state conditions are attained 
where there are no aging effects. In the following sections, 
we analyze the same data for dynamical heterogeneities. 



III. HETEROGENEITY ANALYSIS 

Spatiotemporally heterogeneous dynamics may be 
characterized by a four-point dynamic susceptibility, 
X4(r), which exhibits a peak that grows in proportion 
to the size of the heterogeneities. See for example 
Refs. O [6] for reviews. This method begins with an 
ensemble- averaged self-overlap order parameter, Qt(r), 
constructed so that the contribution from each particle 
decays from 1 toward as time increases from t to t -\- r 
and the particle moves some prescribed distance. If all 
particles experience the same dynamics, then the decay of 
Qti^) vs r will be independent of t. But if the dynamics 
are heterogeneous, then the decay of Qti^) will be faster 
or slower than the time average Q{r) = {Qt{r)) accord- 
ing to the number of fast mobile regions that happen to 
exist at a particular instant. Since the mobile regions are 
independent, this is governed by counting statistics and 



the variance 

X,{T)^N[{QUr))-{Qt{T)f] (1) 

is independent of the number N of particles in the sys- 
tem. The average number n* of particles in a fast mobile 
region has been explicitly computed in Ref. [8 as 



{Qi - QoKQi - Q*) 

where xl is the peak height X4:{'t^)^ Q* is the value of 
Q{t*)^ and Qi and Qo are respectively the average val- 
ues for the slow (Q > Q*) and fast (Q < Q*) regions. 
The same results for n* were found for three very differ- 
ent choices for overlap order parameters, whose associ- 
ated susceptibilities had different peak heights and peak 
times: step function, persistent area, persistent bond [8]. 
Therefore, as long as the prescription of Eq. ([2| is fol- 
lowed, the choice of overlap order parameter is not cru- 
cial. 

Since the video data (eg as in Fig. |2| have insufficient 
resolution to track individual particle positions, we adopt 
an overlap order parameter similar to that introduced in 
Ref. [18 based on image correlations. In particular, we 
divide the video images into 50 narrow strips of constant 
speed and strain rate, 7 pixels ^0.5 jiia wide, containing 
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^ 200 particles. For each of these strips we compute 



Qt(T) 



{iim+,i{t+T))-{my 



(3) 



where (• • •) is the ensemble average over all pixels i run- 
ning along the strip, where di = vr/l, and where / is 
the pixel size. Prior to this, the speed v of the strip 
was found by varying di at fixed r, and averaging over 
t, to maximize the cross-correlation as in the usual PIV 
method. Note that the length scale probed by the asso- 
ciated four-point susceptibility is set by the particle-size 
grayscale variations in the video images. Therefore the 
time r* at which X4.{^) reaches its peak is a characteris- 
tic relaxation time needed for particle-scale relative dis- 
placements. For illustration, example results for Qt{T) vs 
r are shown in Fig. [3^ for a strip corresponding to pack- 
ing fraction ^ = 0.643 and strain rate 7 = 0.0025 s~^. 
This is close to jamming, and indeed the decay is quite 
variable. Multiplying the variance by N gives the suscep- 
tibility shown in Fig. ^jp. This exhibits a peak at delay 
time r* ^ 1000 s of height X4 ^ 9, when the average 
overlap order parameter is Q* ^ 0.5. For a second ex- 
ample strip with a higher strain rate, 7 = 0.078 s~^, the 
dynamics are more homogeneous as seen in Fig. [3] by the 
tighter spread of Qt(r) and the smaller susceptibility. 

To deduce the number n* of particles in a fast mobile 
region from the peak height xl using Eq. ([2|, we must 
first find the three different averages of the overlap order 
parameter. These are shown in Fig.[4]as a function of vol- 
ume fraction, where each point corresponds to a different 
strip and hence to a different strain rate. Since there is no 
evident variation with volume fraction, or strain rate, we 
simply compute a total average over all conditions. The 
average overlap order parameter at the time r* when 
X4 peaks is found to be Q* = Q(r*) = 0.49 ± 0.06. 
The average for the slow regions, where Qt(r*) > Q*, 
is Qi = 0.73 ± 0.08. The average for the fast regions, 
where Qt(r*) < Q*, is Qo = 0.26 ±0.05. These three av- 
erages are indicated by dashed horizontal lines in Fig. |4] 

Since the overlap order parameter is measured in a long 
strip, about half a particle wide and 200 particles long, 
the field of view may not contain the entirety of any of the 
fast moving regions. However if the heterogeneities are 1- 
dimensional string-like swirls, as expected for a quiescent 
system, then they will cut the field of view a number 
of times in proportion to their length. Therefore, the 
true number of particles involved would be a constant 
multiplicative factor larger than the n* value we deduce 
by the above prescription. This is borne out by the good 
comparison of the size of heterogeneities in a monolayer 
of air-fluidized beads, analyzed across the whole sample 
or in strips [32] . But it could also be that in channel flow 
the heterogenieties are linear chains or sheets of particles 
aligned with the velocity, in which the true number of 
particles involved would still be proportional to the n* 
value we deduce. 
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FIG. 5: (Color online) (a) Relaxation time and (b) number of 
particles in a fast-moving heterogeneity, plotted versus strain 
rate, for several volume fractions (f) as labelled. The solid lines 
are fits to a power of 7, with exponents of —1 in (a) and of 
-0.3 in (b). 



IV. RESULTS 

Data for the relaxation time r* and the number n* 
of particles in a fast-moving region are plotted vs strain 
rate in Figs. |5^-b, respectively, for the ~ 1 /im diameter 
colloidal microgel particles. There, each data set corre- 
sponds to a given packing fraction as labelled, and each 
data point corresponds to a different strip and hence to 
a different strain rate. As the strain rate decreases and 
jamming is approached, both r* and n* grow as pow- 
ers of the strain rate. For the time scale, the power 
law is consistent with r* oc I/7 as shown by the solid 
lines. This is the simplest dimensionally-correct possibil- 
ity. For the size scale, all power law fits are consistent 
with n* (X i/^0-30±o.03 ^j^^g exponent agrees with recent 
observations of a value near 1/3 in dry granular systems, 
including experiments on heap flow [18 , simulations of 
uniform shear [33 , and simulations of flow down an in- 
cline [34]. In these works, shear occurs at essentially fixed 
packing fraction near A value in the range 0.2 — 0.3 
was also reported for a Lennard- Jones system [35]. So 
our observations considerably reinforce and extend all 
these results, not just to overdamped systems but also 
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FIG. 6: (Color online) Coefficients of the fits in Fig. [sjto (a) 
T* oc 1/7 and (b) n* oc 1/7°"^ vs volume fraction, shown as 
closed symbols, for the ^ 1 /im diameter particles. The open 
symbols are for the smaller ^0.6 /xm diameter particles. The 
solid curves are fits to a power of Acj) = 10 — (/)c|, where (/)c = 
0.635 is the random close packing fraction, with exponents of 
-4 in (a) and of -1.2 in (b). 
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FIG. 7: (Color online) Scaling collapse of (a) dimensionless 
relaxation time and (b) number of particles in a fast-moving 
heterogeneity versus dimensionless strain rate times Acj) = 
10 — (/)c| to the fourth. Here 770 is the viscosity of water and E 
is the Young elastic modulus of the particulate material. The 
symbol types are the same as in Fig. [5] 



to packing fractions away from 0c both above and below. 
Nonetheless, the value of approximately 1/3 has yet to 
be explained. 

The relaxation time and number of particles in a fast- 
moving region also depend on packing fraction, as well as 
strain rate. This can be seen already in Figs.[5^-b, where 
the data sets shift up and then down as goes from below 
to above 0c- This is displayed more clearly in Fig. [6^-b, 
where the coefficients r*7 and n*7^"^ of the power law fits 
in the previous figure are plotted vs 0. The results grow 
without apparent bound as 0c is approached from either 
side. Data for the smaller particles are also included, 
and display the same behavior. These divergences are 
well described by fits to power laws in A0 = |0 — 0c|, as 
shown, giving r* ex l/A^"^-^^^-^ and n* ex l/A0i-2±O-4. 

Altogether we thus find that the relaxation time and 
the size of fast-moving heterogeneities grow as the jam- 
ming transition is approached as functions of strain rate 
and packing fraction as 

r* cc (7A0^)-\ (4) 
n* oc (7A0^)-°-^ (5) 

Notice that the combination (7A0^) controls the behav- 
ior in both cases. Hence there is more sensitivity to vari- 
ation of A0 than to variation of strain rate. This is qual- 
itatively consistent with numerical results for a driven 



kinetically constrained jamming model [11 . To empha- 
size this feature, we plot all r* and n* results versus 
(7A0^) in Figs. [7^-b, where 7 is rendered dimension- 
less by the intrinsic time scale set by the ratio r]o/E of 
liquid viscosity to particle modulus. Note that this col- 
lapses the data onto power laws with exponents —1 and 
—0.3, respectively. Thus, there are only three exponents 
to explain rather than four. As discussed already, the 
— 1 makes dimensional sense and the —0.3 extends prior 
observations but is not understood. The remaining ex- 
ponent, 4, is reminiscent of the exponent F = 4 in the 
timescale r]o/{EA(j)^) used in Olsson-Teitel scaling plots 
of the shear rheology [29] . 

We note that some simulation results are inconsis- 
tent with Eqs. (|4][5|. In Ref. [3F a Lennard- Jones sys- 
tems shows the same size scaling, xl ^ 1/7^"^ 7 but 
the relaxation time scales quite differently, r* ~ 1/7^-^. 
In Ref. [36] harmonically-repulsive particles under qua- 
sistatic shear show xl 1/A0^-^; the system at nonzero 
strain rates [37] shows xl ^ ]^/^o -5-0.7^ though the de- 
pendence on strain rate may not be a power law. 
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FIG. 8: (Color online) Number n* of particles in a fast-moving 
heterogeneity vs relaxation time r*, made dimensionless by 
the Young modulus E of the particulate material and the 
viscosity 770 of the suspending fluid (water) . The dashed black 
line is a power-law fit with exponent A = 0.31=b0.03as labeled. 
The symbol types are the same as in Fig. [5] 



V. DISCUSSION 

The observations made here, summarized by Eqs. ^ 
[5|, combine to give the size of heterogeneities as a power- 
law of the relaxation time: 



* \ A 



n* oc (r*) 



(6) 



For emphasis, we plot n* data vs r^E/rfo in Fig. [s] on 
logarithmic axes and observe that indeed the data col- 
lapse to a straight line. Fitting to a power-law gives the 
exponent and uncertainty as 



A = 0.31 ±0.03. 



(7) 



Note that the dynamic range is more than two decades 
in size and seven in relaxation time, which is sufficient to 
rule out the possibilities of a logarithmic or exponential 
connection between n* and r*. The power law form is 
consistent with mode-coupling theory, and the exponent 
is only slightly smaller than the expectation I/7 reported 
in Refs. [9 , 10 . 

One advantage of plotting n* and r* parametrically 
versus one another, rather than versus the control pa- 
rameters, is that it allows comparison with other systems 
where there is no shear or where the control parame- 
ter is something other than strain rate. For example, 
A = 1 /2 is reported for the driven kinetically constrained 
model mentioned above [11 , while a logarithmic con- 
nection better accounts for the simulations of Brownian 
harmonically-repulsive particles [38]. In terms of experi- 
ment, comparison is possible for only a few experiments 
of which we are aware. For hard spheres, Ref. [12 found 
that r* grows faster than a power law, and that x% grows 
slower than a power law, as (j) approaches (j)c from below. 
It is stated that n* grows logarithmically with r*, which 
disagrees with our results. For soft NIPA microgel par- 
ticles similar to those studied here, Ref. [27 found that 



both r* and xl grow with increasing packing fraction; 
no functional form was proposed or tested. We digitized 
their data and plot parametrically, rather than vs pack- 
ing fraction, and find power law behavior of the form 
Eq. (|6| with exponent A = 0.34 ± 0.16. This is consis- 
tent with our findings. For a monolayer of large spher- 
ical grains fluidized by a steady upflow of air, Ref. [17] 
found that there is a meaningful effective temperature 
Teff, and that size and time scales are consistent with 
n* - l/Teff°-^^°-^ and r* - l/Teff^^°•^ respectively. 
These combine to give a power law relationship, Eq. (|6|, 
with exponent A = 0.35 ±0.15 that agrees with the re- 
sults here in Fig. [8] Further experiments on the fluidized 
grains, where the sample is tilted and where the analy- 
sis is carried out by dividing the sample into a series of 
strips each at a different pressure, also appear to agree 
[32] . For steady gravity-driven flow of grains down along 
a confined heap, and visualized through the sidewalls as a 
function of depth z below the free surface, Ref. [18 found 
r* (X l/I and n* cx (l//)'^ where / = '^dj -^Jgz is the in- 
ertia number, d is the grain diameter, g = 9.8 m/s^, and 
A = 0.33 ±0.02. This system is underdamped, so the time 
scale is rendered dimensionless by different microscopic 
physics, but the exponent for the power-law connection 
between n* and r* is the same as found here. 



VI. CONCLUSION 

In this paper we presented a study of dynamical het- 
erogeneities in a colloidal system that (a) is not hard 
spheres, (b) is compressed above as well as below and 
(c) is subject to shear. Our experiments made crucial 
use of custom synthesized NIPA microgel particles and 
of a custom fabricated microfluidic channel, as well as 
of a novel video-based dynamical order parameter. As 
jamming is approached by bringing the packing fraction 
difference A<p = |^ — or the strain rate 7, to zero, we 
demonstrated that the time and size scales for dynamical 
heterogeneities both grow aspowers of the combination 
(7A(/)^) according to Eqs. (|4|5|. While there is prece- 
dence for the observed strain rate dependence from ex- 
periments on underdamped granular systems, the pack- 
ing fraction dependence appears to be a new result. The 
observed connection between the size and time scale is 
a power law, n* oc (r*)'^, consistent with mode coupling 
theories but not with observations for an unsheared sus- 
pension of Brownian hard sphere colloids. It is intriguing 
that the exponent we find, A ~ 1/3, agrees with prece- 
dents for unsheared soft particles [27] and for two driven 
systems of hard grains - one with shear [T8^ and one with- 
out [TTl [32] . This suggests universality with respect to 
interactions, but with unsheared hard sphere colloids in 
a different universality class. 
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